#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Difference Map                                                     #
#                                                                           #
# (C) 2dx.org, GNU Plublic License.                                         #
#                                                                           #
# Created..........: 08/02/2012                                             #
# Last Modification: 08/02/2012                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 70 
#
# MANUAL: This script makes use of CCP4 commands to generate difference map.
#
# REMARK: The following section defines variables that will be displayed in the GUI:
# DISPLAY: RESMAX
# DISPLAY: RESMIN
# DISPLAY: realcell 
# DISPLAY: realang
# DISPLAY: diff_map1
# DISPLAY: diff_map1_name
# DISPLAY: diff_map2
# DISPLAY: diff_map2_name
# DISPLAY: diffmap_cntrs_mode
# DISPLAY: diffmap_cntrs_max
# DISPLAY: diffmap_cntrs_min
# DISPLAY: diffmap_cntrs_step
# DISPLAY: npo_cntrs_step 
# DISPLAY: npo_cntrs_below 
#
# REMARK: The following section defines the variables that the GUI will fill in from the database, before launching this script:
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set imagename = ""
set tempkeep = ""
set RESMIN = ""
set RESMAX = ""
set realang = ""
set realcell = ""
set SYM = ""
set ctf_ttf = ""
set diff_map1 = ""
set diff_map1_name = ""
set diff_map2 = ""
set diff_map2_name = ""
set diffmap_cntrs_mode = ""
set diffmap_cntrs_max = ""
set diffmap_cntrs_min = ""
set diffmap_cntrs_step = ""
set npo_cntrs_below = ""
set npo_cntrs_step = ""
#
#$end_vars
#
echo "<<@progress: 1>>"
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
set scriptname = 2dx_diffmap
#
set IS_2DX = yes
source ${proc_2dx}/2dx_makedirs 
#
\rm -f LOGS/${scriptname}.results
#
echo "<<@evaluate>>"
#
set date = `date`
echo date = ${date}
#
set filename1 = `echo $diff_map1 | awk -F/ '{print $NF }' | awk -F. '{ print $1 }'`
set filename1 = "diff_map1_${filename1}"
echo "the filename of map 1: ${filename1}"
set filename2 = `echo $diff_map2 | awk -F/ '{print $NF }' | awk -F. '{ print $1 }'`
set filename2 = "diff_map2_${filename2}"
echo "the filename of map 2: ${filename2}"
#
##############################################################################
${proc_2dx}/linblock  "Scaling amplitudes" 
##############################################################################
if ( -e ${diff_map1} ) then
	echo "scaling amplitudes of  ${diff_map1} "
	\rm -f  ${filename1}_scaled_MRClefthanded.mtz  
	${proc_2dx}/2dx_scaleamp_ccp4.com ${diff_map1} ${filename1}_scaled_MRClefthanded.mtz
	echo "# IMAGE: ${filename1}_scaled_MRClefthanded.mtz <MTZ: Map 1 with scaled amplitudes >" >> LOGS/${scriptname}.results
else
	${proc_2dx}/protest "Could not find ${diff_map2}"
endif
if ( -e ${diff_map2} ) then
	echo "scaling amplitudes of  ${diff_map2} "
	\rm -f  ${filename2}_scaled_MRClefthanded.mtz  
	${proc_2dx}/2dx_scaleamp_ccp4.com ${diff_map2} ${filename2}_scaled_MRClefthanded.mtz
	echo "# IMAGE: ${filename2}_scaled_MRClefthanded.mtz <MTZ: Map 2 with scaled amplitudes >" >> LOGS/${scriptname}.results
else
	${proc_2dx}/protest "Could not find ${diff_map2}"
endif
echo "<<@progress: 10>>"
##############################################################################
${proc_2dx}/linblock  "Limiting the resolution to ${RESMAX}" 
##############################################################################
\rm -f  ${filename1}_res_MRClefthanded.mtz  
${bin_ccp4}/sftools << eot
read ${filename1}_scaled_MRClefthanded.mtz
select resolution >= ${RESMAX}
checkhkl
write ${filename1}_res_MRClefthanded.mtz
end
eot
#
echo "# IMAGE: ${filename1}_res_MRClefthanded.mtz <MTZ: Map 1 (resolution limited) >" >> LOGS/${scriptname}.results
#
\rm -f  ${filename2}_res_MRClefthanded.mtz  
${bin_ccp4}/sftools << eot
read ${filename2}_scaled_MRClefthanded.mtz
select resolution >= ${RESMAX}
checkhkl
write ${filename2}_res_MRClefthanded.mtz
end
eot
#
echo "# IMAGE: ${filename2}_res_MRClefthanded.mtz <MTZ: Map 2 (resolution limited) >" >> LOGS/${scriptname}.results
#
echo "<<@progress: 20>>"
#
##############################################################################
${proc_2dx}/linblock  "Inverse fft to get map" 
##############################################################################
set realcell_a = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = int($1)  } END { print s }'`
set realcell_b = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = int($2)  } END { print s }'`
set realang2 = `echo ${realang} | sed 's/,/ /g' | awk '{ s = int($1)*2  } END { print s }'`
set grid_a = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = (int($1 -  ($1 % 4)))*2  } END { print s }'`
set grid_b = `echo ${realcell} | sed 's/,/ /g' | awk '{ s = (int($2 -  ($2 % 4)))*2  } END { print s }'`

if( $grid_a  < $realcell_a || $grid_b < $realcell_b ) then
	echo  ":: The unit cell was cropped to ${grid_a}, ${grid_b}, ${realang}. (Dividable by 4)" 
endif
#
\rm -f ${filename1}.map
#
${bin_ccp4}/fft HKLIN   ${filename1}_res_MRClefthanded.mtz MAPOUT  ${filename1}.map \
    << eot
resol 80 ${RESMAX}
! Fo
TITLE "${RESMAX} EM map for ${filename1}"
LABIN F1=F PHI=PHI W=FOM
scale F1 1.0
AXIS Y X Z
GRID ${grid_a} ${grid_b} ${realang2}
end
eot
#
\rm -f ${filename2}.map
#
${bin_ccp4}/fft HKLIN   ${filename2}_res_MRClefthanded.mtz MAPOUT  ${filename2}.map \
    << eot
resol 200 ${RESMAX}
! Fo
TITLE "${RESMAX} EM map for ${filename2}"
LABIN F1=F PHI=PHI W=FOM
scale F1 1.0
AXIS Y X Z
GRID ${grid_a} ${grid_b} ${realang2} 
end
eot
#
echo "<<@progress: 30>>"
##############################################################################
${proc_2dx}/linblock  "npo - to create a line plot for both maps" 
##############################################################################
${bin_ccp4}/npo  MAPIN ${filename1}.map PLOT ${filename1}.plot <<EOF
TITLE   ${diff_map1_name}
MAP SCALE 0.5
CONTRS SIG -3.0  to 3.0 BY ${npo_cntrs_step} 
MODE BELOW 0 DASHED 1 0.25 0
MODE BELOW ${npo_cntrs_below}  DASHED 1 0.25 0
SECTS 0 0
PLOT 
EOF
#
${bin_ccp4}/npo  MAPIN ${filename2}.map PLOT ${filename2}.plot <<EOF
TITLE   ${diff_map2_name}
MAP SCALE 0.5
CONTRS SIG -3.0  to 3.0 BY ${npo_cntrs_step} 
MODE BELOW 0 DASHED 1 0.25 0
MODE BELOW ${npo_cntrs_below}  DASHED 1 0.25 0
SECTS 0 0
PLOT 
EOF
#
echo "<<@progress: 35>>"
#############################################################################
${proc_2dx}/linblock "laserplot - to create PS/${diff_map1_name}.ps"
#############################################################################
#
\rm -f PS/${diff_map1_name}.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${diff_map1_name}.ps ${filename1}.plot
\rm -f ${filename1}.plot
echo "# IMAGE-IMPORTANT: PS/${diff_map1_name}.ps <PS: map 1>"  >> LOGS/${scriptname}.results
#
\rm -f PS/${diff_map2_name}.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${diff_map2_name}.ps ${filename2}.plot
\rm -f ${filename2}.plot
echo "# IMAGE-IMPORTANT: PS/${diff_map2_name}.ps <PS: map 2>"  >> LOGS/${scriptname}.results
#
echo "<<@progress: 40>>"
##############################################################################
${proc_2dx}/linblock  "npo - to plot the overlay" 
##############################################################################
${bin_ccp4}/npo  MAPIN ${filename1}.map MAPIN2 ${filename2}.map PLOT overlay.plot <<EOF
TITLE ${diff_map1_name} & ${diff_map2_name} overlay
MAP SCALE 0.5 
CONTRS SIG -2.0 to 2.0 BY ${npo_cntrs_step}	  #1st map  
MODE BELOW  ${npo_cntrs_below} DASHED 1 0.25 0
MODE RED
CONTRS SIG -2.0 to 2.0 BY ${npo_cntrs_step}       #2nd map
MODE GREEN
MODE BELOW  ${npo_cntrs_below} DASHED 1 0.25 0
SECTS 0 0
PLOT 
EOF
#
echo "<<@progress: 45>>"
#############################################################################
${proc_2dx}/linblock "laserplot - to create PS/${diff_map1_name}_${diff_map2_name}_overlay.ps"
#############################################################################
\rm -f PS/${diff_map1_name}_${diff_map2_name}_overlay.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${diff_map1_name}_${diff_map2_name}_overlay.ps overlay.plot
\rm -f overlay.plot
echo "# IMAGE-IMPORTANT: PS/${diff_map1_name}_${diff_map2_name}_overlay.ps <PS: overlay of both maps >"  >> LOGS/${scriptname}.results
#
echo "<<@progress: 50>>"
#############################################################################
${proc_2dx}/linblock "overlapmap -  to calulate difference map - ${diff_map1_name} - ${diff_map2_name}"
#############################################################################
${bin_ccp4}/overlapmap mapin1 ${filename1}.map mapin2 ${filename2}.map mapout diffmap_${filename1}-${filename2}.map <<eof
map add 1 -1
end 
eof
#
#
echo "<<@progress: 60>>"
##############################################################################
${proc_2dx}/linblock  "npo - to create a line plot for difference map" 
##############################################################################
${bin_ccp4}/npo  MAPIN diffmap_${filename1}-${filename2}.map  PLOT diffmap.plot <<EOF
TITLE Difference Map of ${diff_map1_name} - ${diff_map2_name} 
MAP SCALE 0.5
CONTRS ${diffmap_cntrs_mode} ${diffmap_cntrs_min} to ${diffmap_cntrs_max}  BY ${diffmap_cntrs_step}
MODE GREEN #BLUE 		#diff_map1 
MODE BELOW 0 YELLOW #RED 	#diff_map2 
SECTS 0 0
PLOT 
EOF
#
#
echo "<<@progress: 70>>"
#############################################################################
${proc_2dx}/linblock "laserplot - to create  PS/${filename1}-${filename2}_diffmap.ps"
#############################################################################
\rm -f PS/${filename1}-${filename2}_diffmap.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${filename1}-${filename2}_diffmap.ps diffmap.plot
\rm -f diffmap.plot
echo "# IMAGE-IMPORTANT: PS/${filename1}-${filename2}_diffmap.ps <PS: Difference Map >"  >> LOGS/${scriptname}.results
#
echo "<<@progress: 80>>"
##############################################################################
${proc_2dx}/linblock  "npo - to create a line plot for map 1 & difference map" 
##############################################################################
${bin_ccp4}/npo  MAPIN ${filename1}.map MAPIN2 diffmap_${filename1}-${filename2}.map PLOT diff_map1_diff.plot <<EOF
TITLE ${diff_map1_name} and the difference map
MAP SCALE 0.5
CONTRS SIG -3.0  to 3.0 BY ${npo_cntrs_step}	#1st map 
MODE BELOW 0 DASHED 1 0.25 0
CONTRS ${diffmap_cntrs_mode} ${diffmap_cntrs_min} to ${diffmap_cntrs_max}  BY ${diffmap_cntrs_step}	#2nd map 
MODE YELLOW
MODE BELOW 0 GREEN
SECTS 0 0
PLOT 
EOF
#
echo "<<@progress: 90>>"
#############################################################################
${proc_2dx}/linblock "laserplot - to create PS/${filename1}and${filename1}-${filename2}_diffmap.ps"
#############################################################################
\rm -f PS/${filename1}and${filename1}-${filename2}_diffmap.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${filename1}and${filename1}-${filename2}_diffmap.ps diff_map1_diff.plot
\rm -f diff_map1_diff.plot
echo "# IMAGE-IMPORTANT: PS/${filename1}and${filename1}-${filename2}_diffmap.ps <PS: <Map 1 and Difference Map >"  >> LOGS/${scriptname}.results
#
###########################################################################
#############################################################################
#
echo "<<@progress: 100>>"
#
##########################################################################
${proc_2dx}/linblock "${scriptname} - normal end."
##########################################################################
#
exit
#
## Just to make sure it is displayed in the 2dx GUI:
#source ${proc_2dx}/2dx_checklattice_sub.com 
##
##
